A message passing approach for general epidemic models 
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In most models of the spread of disease over contact networks it is assumed that the probabilities 
per unit time of disease transmission and recovery from disease are constant, implying exponential 
distributions of the time intervals for transmission and recovery. Time intervals for real diseases, 
however, have distributions that in most cases are far from exponential, which leads to disagree- 
ments, both qualitative and quantitative, with the models. In this paper, we study a generalized 
version of the SIR (susceptible-infected-recovered) model of epidemic disease that allows for arbi- 
trary distributions of transmission and recovery times. Standard differential equation approaches 
cannot be used for this generalized model, but we show that the problem can be reformulated as 
a time-dependent message passing calculation on the appropriate contact network. The calculation 
is exact on trees (i.e., loopless networks) or locally tree-like networks (such as random graphs) in 
the large system size limit. On non-tree-like networks we show that the calculation gives a rigor- 
ous bound on the size of disease outbreaks. We demonstrate the method with applications to two 
specific models and the results compare favorably with numerical simulations. 



I. INTRODUCTION 

The mathematical modeling of infectious disease out- 
breaks in human populations has a long history, stretch- 
ing back to the pioneering work of Lowell Reed, Ander- 
son McKendrick, and others in the early twentieth cen- 
tury [1] . The standard analytic approach involves divid- 
ing the modeled population into classes or compartments 
according to their status with respect to the disease of 
interest — uninfected but susceptible, infected, recovered, 
and so forth — and then writing differential equations to 
describe the mass flow of individuals between compart- 
ments according to the dynamics of the infection pro- 
cess [1, 2], 

Such compartmental models have proven flexible, 
tractable, and highly informative as a general guide to the 
population-level behavior of diseases, but they also suffer 
from a number of serious deficiencies, of which two are 
particularly significant. The first, which has attracted a 
lot of recent attention in the literature, is the assumption 
of random mixing. In order to write differential equations 
for flows between compartments, we must make a fully 
mixed or mass-action approximation whereby we assume 
that the probability of disease-causing contact with any 
member of a particular compartment is the same. In real 
life this is far from true — most people have high prob- 
ability of contact with only that small fraction of the 
population they rub shoulders with regularly, and a very 
small chance of contact with everyone else. The incorpo- 
ration of more realistic mixing patterns into epidemiolog- 
ical modeling has given rise to the field of network epi- 
demiology, in which contacts are modeled as a network, 
either static [3-9] or dynamic [10-12], and the structure 
of the network can have a profound impact on the spread 
of the disease [13-16]. 

In this paper, however, we focus on a different short- 
coming of compartmental models, one that has by com- 
parison received little attention, but which is at least as 
important as the mass-action approximation. In order to 



write down the differential equations of a compartmental 
disease model, one must make the assumption that move- 
ment between compartments takes place at a stochasti- 
cally constant rate. In modeling a disease from which 
most victims recover, for instance, one typically assumes 
that an infected individual has a constant probability per 
unit time of recovery. While being a useful assumption 
from a mathematical point of view, however, this be- 
havior is very far from that of most real diseases. The 
assumption of constant probability of recovery implies an 
exponential distribution of times for which individuals re- 
main infected, so that the most probable duration of in- 
fection is zero, and probability decreases uniformly with 
time. In reality, most diseases show a roughly constant 
duration of infection — a week, say, or a month — with rel- 
atively small fluctuations from person to person, so that 
the distribution of durations has a sharp peak about the 
average value and is highly nonexponential. Such nonex- 
ponential distributions are known to have a substantial 
effect, both qualitative and quantitative, on the shape of 
epidemics [17-21]. 

If one is willing to make the mass-action approxi- 
mation, then nonexponential behavior can be incorpo- 
rated into epidemic models by reformulating the the- 
ory in terms of integro-differential equations [22, 23]. If, 
however, one wishes also to retain the advances of net- 
work epidemiology in representing nonrandom contact 
patterns, then even this approach does not work and a 
new method of solution is necessary. In this paper we 
demonstrate that in the latter case the calculations can 
be reformulated in the language of message passing algo- 
rithms of the kind known as belief propagation or sum- 
product methods. In addition to providing exact solu- 
tions for the dynamics of quite general epidemic models 
on large classes of networks, the message passing formu- 
lation also leads to a number of other results concerning 
network epidemiology, including a rigorous upper bound 
on the size of disease outbreaks, results for late-time be- 
havior, and results for the average behavior of epidemics 
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in random network ensembles. 
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II. A MESSAGE PASSING FORMULATION OF 
EPIDEMICS 

We begin by denning the problem. We consider the 
simplest nontrivial model of epidemic disease, the SIR 
model, in which an individual can be in one of three 
disease states, susceptible, infected, or recovered. We 
will assume an initial condition for the epidemic in which 
each vertex is susceptible with independent probability z 
and infected otherwise. 

We assume that disease transmission is taking place on 
a given contact network, meaning that disease can only 
be transmitted between individuals who are directly con- 
nected by an edge in that network. We also generalize 
the model to allow for nonexponential distributions of 
the times at which transitions between these states oc- 
cur, i.e., the times at which infection and recovery occur. 
To be completely general, let us define s(r) dr to be the 
probability that an individual infected with the disease of 
interest first makes contact sufficient to transmit the dis- 
ease to a particular network neighbor at a time between 
r and t + dr after their infection. Similarly let us define 
r(r) dr to be the probability that an individual infected 
with the disease recovers from it at a time between r and 
t + dr after infection. 

An infected individual can only transmit the disease to 
a susceptible neighbor if they are still infected at the time 
of contact, and hence the probability that transmission 
actually occurs between t and r + dr after infection is 
equal to the probability s(r) dr times the probability 
J T °° t{t') dr' that the individual has not yet recovered. 
Let us denote this overall probability of transmission by 
f(r) dr: 



/(r) dr = s(t) dr / r(r') dr'. 



(1) 



Note that this function, unlike s(r) and r(r) does not in- 
tegrate to unity. Rather, it integrates to the total proba- 
bility that a vertex transmits the disease to its neighbor 
before it recovers, a probability referred to elsewhere var- 
iously as the transmissibility or infectivity of the disease. 

The fundamental quantity appearing in our message 
passing formulation of disease transmission — the "mes- 
sage" that is passed among network vertices in the 
calculation — is the probability, which we denote H l ^i (t) , 
that a vertex j has not passed the disease to neighboring 
vertex i by absolute time t. (Without loss of general- 
ity, we will assume the epidemic to begin at absolute 
time t = 0.) An especially simple case of our approach 
arises when the network of interest takes the form of a 
tree, i.e., a network having no loops. In this case, the 
failure of vertex j to pass the disease to vertex i before 
time t can occur in either of two ways, as illustrated in 
Fig. 1. First, it may be that, if and when vertex j con- 
tracts the disease, it fails to transmit it to i within an 
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FIG. 1: The probability that vertex i does not contract the 
disease from its neighbor j before time t is equal to the prob- 
ability that j fails to transmit the disease within an interval t 
of catching it, plus the probability that it does transmit the 
disease within an interval t but that j received the disease 
from its neighbors (here denoted k, q, and r) too late to pass 
it on to i in time. 



interval t from infection, in which case clearly i does not 
contract the disease before absolute time t. The proba- 
bility of this occurrence is 1 — f Q /(r) dr. 

The second possibility is that j is scheduled to trans- 
mit the disease within time t of contracting it, but that 
j itself got the disease (from one of its other neighbors) 
too late for that transmission to occur before absolute 
time t, or indeed never got the disease at all. If j trans- 
mits the disease at time r after contracting it, but fails 
to contract the disease before time t' = t — t then i 
does not receive the disease before time t. The prob- 
ability that j does not contract the disease before t' is 
z YlieAf(j)\i H^ l (t') 7 where the leading factor of z repre- 
sents the probability that j was not one of those vertices 
initially infected with the disease at t = 0. The notation 
J\f(j)\i denotes the set of neighbors of vertex j, exclud- 
ing vertex i. Now integrating over t', we find the total 
probability that j fails to transmit the disease before t to 
be zf*f(t-t>)ll lemj)y H^ l (t')dt'. 

Putting the two contributions to fP <_J (t) together and 
writing t — t' = r we arrive at the message passing equa- 
tion 



H^(t) = l- ff(r)\l-z J] H^\t-r) 



dr. (2) 



For the special case of a network taking the form of a 
tree, this equation gives us, at least in principle, a com- 
plete solution for the probabilities H l< ^i{t) for all t and 
arbitrary /(r). 

Normally, however, H l< ~ J is not the quantity of epi- 
demiological interest. More commonly one wants to know 
things such as the fraction of the population that will be 
in the various disease states at different times, or more 
generally the probability that a particular individual will 
be in each disease state. Let us define Si(t) to be 1 if 
individual i is susceptible at time t and otherwise, and 
similarly define h (t) and Ri (t) for the infected and recov- 
ered states. Then P(Si(t) = 1) denotes the probability 
that vertex i is susceptible at time t. For the sake of 
economy we will also write this probability more briefly 
simply as P{Si). For i to be susceptible at time t we re- 
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quire (a) that i is not one of the vertices initially infected 
at t = 0, which happens with probability z, and (b) that 
i not receive the infection from any of its neighbors be- 
fore time t. Thus P(Si) can be expressed quite simply 



as 



P(S l ) = z 1 [[ H^(t). 



(3) 



Once we have P(<Sj), however, one can also imme- 
diately calculate P(Ii) and P(Ri). Note that the rate 
dP(Ii)/dt at which P{Ii) increases is equal to the rate 
at which P(Si) decreases — since all individuals moving 
out of the susceptible state must move into the infected 
state — minus the rate at which i recovers. The recovery 
rate has two contributions: the probability 1 — z that i 
was infected at time t = times the rate r(t) of recovery 
a time t later, and the probability that i was infected 
at some later time t' < t (which is simply — dP{Si)/dt') 
times the rate r(t — t') of recovery a time t — t' later. This 
allows us to write a rate equation for P{h) thus: 



dP(Ii) dP{Si) 



dt 



dt 



-(l-z)r(t)- 



1 , dP(SA , . 



(4) 

By integrating this equation we can calculate P{h) for 
any t, and then we can calculate P(Ri) from the knowl- 
edge that the probabilities of the three states must sum 
to one: 



P(Ri) = 1 - P(Si) - P{h). 



(5) 



Between them, Eqs. (2) to (5) now give us a complete so- 
lution for the three probabilities, for arbitrary (including 
nonexponential) distributions of infection and recovery 
times. 



III. MESSAGE PASSING ON NON-TREE 
NETWORKS 

The developments of the previous section give us a so- 
lution for the SIR model in the case where the network 
of interest has no loops, but almost all real-world net- 
works have loops, and usually many of them. It is known 
that message passing methods, while not exact on non- 
tree networks can still give good approximate answers 
in some cases. In the present case, however, we can go 
further than such qualitative statements and show that 
our message passing calculation provides a rigorous upper 
bound to the number of infected individuals on networks 
that contain loops. To prove this result, consider the fol- 
lowing alternative formulation of the epidemic process. 

In the generalized SIR model discussed here, evolution 
of the disease involves infected individuals spreading in- 
fection to their susceptible neighbors at times after infec- 
tion drawn from the distribution s(t) and recovering at 
times after infection drawn from r(r). There is, however, 
no requirement that we draw these times at the moment 



of infection. We can if we wish draw them ahead of time, 
before executing the steps of the model. That is, we can 
for each vertex i in the network draw a time Tj from 
the distribution r(r) and associate it with that vertex. 
When vertex i becomes infected, we look up the value 
of Tj which tells us the interval of time before i recovers. 
For the edges the situation is only a little more compli- 
cated. We replace each undirected edge in the network 
with two directed edges pointing in opposite directions, 
to represent the act of disease transmission in either di- 
rection between the two relevant vertices. Then for each 



directed edge j 



we draw a time 



from the dis- 



tribution s(w) to represent the time after infection of j 
at which contact is made with i. If this time falls before 
the recovery of j, i.e., if w^ < Tj, then transmission will 
take place if j is ever infected, and will occur an inter- 
val Wij after infection. If, however, j recovers first, i.e., if 
Wij > Tj, then no transmission takes place, which we can, 
if we wish, represent mathematically by setting w^ = oo. 

The end result is a directed "transmission network" in 
which the edges represent possible transmission events 
and the values w^ on the edges represent the time delay 
between arrival of the infection at j (if that ever happens) 
and arrival of the infection at i. 

In terms of this network it is now quite simple to write 
down the probability P(Si) that vertex i is susceptible at 
time t. In order to be susceptible we require (a) that i was 
not infected at time 0, which happens with probability z, 
and (b) that there exists no path from any vertex j to 
vertex i such that vertex j was infected at time and the 
sum of the time delays Wq along the path is less than t. 

An alternative way of thinking about this second con- 
dition is to consider the neighborhood of radius t about 
vertex i, meaning the set of vertices j a distance t or less 
from i, where distance is measured in terms of the sum 
of the values w^ along the path — the shortest weighted 
distance in the language of graph theory. If any of the 
vertices in this neighborhood is infected at time zero 
then vertex i will not be susceptible at time t. Let us 
suppose that there are rn vertices in the neighborhood, 
excluding vertex i itself. Then the probability that i is 
susceptible — for this particular choice of the w^ and Tj — 
is z ,li+1 . We are interested, however, in the probability 



averaged over all values of the w^ and Tj, which is 

P{Si)=z{z n <), 



(6) 



where the angle brackets (. . .} denote the average over 
the ensemble of values of and Tj. 

This equation is correct and exact in all cases. To re- 
late it to our previous message passing approach and un- 
derstand how the calculation proceeds on networks with 
loops, consider the alternative set of vertex counts ny, 
which are the numbers of vertices whose distance to i is t 
or less, but now with the restriction that the penultimate 
vertex along the path to i must be vertex j. For reasons 
that will shortly become clear, we also forbid paths that 
pass through vertex i more than once. That is, there may 
be a path of length t or less that first passes through i to 



4 




FIG. 2: A small directed transmission network in which each 
edge is labeled with its associated transmission delay Wij , ex- 
cept for edges with toy > r, , which are labeled oo. The three 
red vertices denote those within distance 6 of the black vertex 
and the red edges correspond to the weighted shortest paths 
from the red vertices to the black one. If we approximate 
the number of red vertices as in Eq. (7) by the sum of the 
numbers of vertices within distance 6 that are reachable via 
each of the black vertex's immediate neighbors, then we will 
count four vertices instead of three: the top red vertex will be 
counted twice because the blue edge provides a second path 
from this vertex to the black one. 



reach j and then returns to i, but such paths are disal- 
lowed. In practice, a simple way to enforce this constraint 
is to remove from the network all directed edges outgoing 
from vertex i. In this case, vertex i is said to be a cavity 
vertex or in the cavity state. 

We now observe that, as illustrated in Fig. 2, the sum 
of riij over all neighbors j is always at least as great as nf. 



rii< 2J n i 



(7) 



where the inequality becomes an exact equality if the 
network is a tree. (It is in order to ensure this equality 
that we exclude paths that pass through i twice.) Then 

z m > z £ 3eA r(i) ™« an d 

P{S t ) = z(z n >) > z(z^^> n ^) = z ( J[ z n A. (8) 

We now apply a version of the Chebyshev integral in- 
equality, proved in the appendix, that for any set of non- 
negative functions fi(x±, . . . ,x&), . . . , f n (xi, ■ ■ ■ ,%k) that 
are monotone increasing or decreasing in every argument, 
says 



i=l 



))>l[(fi(xi,...,x k )), (9) 



where the average is over any distribution of independent 
variables Xi, . . . ,Xk- Applied to Eq. (8), this inequality 
tells us that 

p(Si)>z n (z n »)=z n H^(t), (io) 



where we have defined 

H^ j {t) = (z n »). 



This quantity is the average probability that at time t 
the infection has not been passed to vertex i via neigh- 
bor j (again excluding cases where the infection passes 
through i twice). It plays the same role as the corre- 
sponding quantity in Eq. (3) for the case of a tree, and 
we can evaluate it in an analogous way. As before we 
split H % ^if) into two parts. The first is the probability 
that, even if j is infected, it does not transmit the dis- 
ease to i within time t of infection. This probability, as 
before, is 1— J Q /(r)dr, where /(t) is defined by Eq. (1). 

The second part is the probability that j is sched- 
uled to transmit the disease within time r < t of con- 
tracting it, but that j itself gets the disease too late 
for the transmission to occur before absolute time t (or 
j never gets the disease at all). For transmission be- 
fore time t vertex j needs to contract the disease before 
t' = t — t and the probability that this does not hap- 
pen is P(Sj(t')\i in cavity), where it is now important 
that i is in the cavity state, so as to disallow paths for 
infection that pass through i itself. Then the probabil- 
ity that j fails to transmit the disease before time t is 

So fit - t')P(Sj{~t')\i in cavity) dt'. 

The probability P(Sj(t')\i in cavity) we can calculate 
from the appropriate analog of Eq. (10) but with both i 
and j in the cavity state, i.e., with their outgoing edges 
deleted. But consider now adding back in all the edges 
leading from i except the one to j. In doing so we only 
add paths to the network and hence potentially increase 
the size of the neighborhood of vertex j but never de- 
crease it. This implies that we only decrease P(Sj(t')), 
so that 

P(Sj(t')\i in cavity) > P(Sj(t')\i ~> 3 deleted) 

>z <*"*) 

leAT(J)\i 

= z n H^\t'), (12) 

where we have used Eq. (10). Combining our two contri- 
butions to H l ^^{t) and writing t — t' = t, we now find 
that 



H 



(*)>1- 



i- [ m\i-z n H^\t- T ) 



dr. (13) 



This result is very similar to the message passing equal- 
ity of Eq. (2) , but it is an inequality, and hence cannot be 
directly employed to calculate properties of the epidemic. 
Let us, however, define a different function F l ^i(t) by 
the equation 



**-'■(«) = 1- ff(r)\l-z [] F*-'(t-r) 
Jo l 



!SAT(i)\i 



dr, (14) 



(11) 



which is an equality and so can be used to calculate F 14 ^ 3 , 
for instance by iteration starting from a suitable ini- 
tial value -Fq <_,7 (£). Suppose we choose as our initial 
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value Fq J (t) — H l ^ 3 {t) for all i, j and t, so that, from 



Eq. (13), we have 
J^ J '(t)>l-/Wl-* II Ft\t-r) 



l£SS(j)\i 



dr. (15) 



(Of course we don't know the value of H 14 ^ 3 (t), but sup- 
pose for the moment that we do.) Then, performing one 
step of the iteration, we arrive at a new value F[^ J (t) 
thus: 



(*) 



1 



f(r) 



1 



- n f o 

ieSS(j)\i 



\t-r) 



<Ft j (t), 



dr 



(16) 



where we have used Eq. (15). But note that, since /(t) > 
for all r, Eq. (16) also implies that 



ff(r) II Fr\t-r)&T 

ieM(j)\i 

</V(r) II Ft l (t-r)dr, (17) 



and hence from Eq. (16) 



(*)>!-/' 
Jo 



f(r) 



i-z n ^r'(i-r) 



dr, (18) 



which is the equivalent of Eq. (15) for F^ 3 (t). Now we 
can repeat the same argument to show that for a general 
step of the iteration we must have 



*r j (*)<-C- j i(*)- 



(19) 



In the limit m — > oo, the iteration must converge, since 
F-£~ l (t) is bounded below by 1 — J* /(r) dr, and hence 
in this limit we get a solution to Eq. (14) that satisfies 

F l< ~ 3 (t) < F^ j (t) = H i4 ~ j (t) , (20) 

for all i,j and t. 

Now, making use of Eq. (10), we have 

P(Si) > z J] H^ 3 (t) > z 11 F^ 3 (t). (21) 

jeAT(i) jeAT(i) 

Thus Eq. (14) allows us to calculate a rigorous lower 
bound on the probability that any vertex is in the sus- 
ceptible state. Notice that Eq. (14) is the same as the 
equation for H 1 ^ 3 in the tree case, Eq. (2), but is per- 
fectly well defined for any network, tree or otherwise. 

Our lower bound on P(Si) also gives us upper bounds 
on P(Ii) and P(Ri), both of which are trivially less than 
1 — P(Si), as well as an upper bound on the sum P{h) + 
P(Ri) = 1 — P(Si), which is the total probability that i 
has ever caught the disease. Hence our message passing 



calculation can in this case give us an upper bound on the 
number of individuals infected by an epidemic, a result 
of possible value — a guarantee that infection will not rise 
above a certain level could be used as a quality function to 
quantify the efficacy of proposed vaccination campaigns 
or other public health interventions. 

Employing Eqs. (14) and (21) in a message passing al- 
gorithm would involve propagating messages that take 
the form of functions F l4 ~ 3 (t) of time. On a tree, one 
would start with the leaves of the tree, for which Eq. (14) 
is trivial, and work inwards through the network until the 
functions on all edges have been evaluated. On a non- 
tree network, the calculation is more complicated because 
one does not in general know any of the F li ~ 3 {t) to be- 
gin with, so one would have to make an initial guess and 
then iterate Eq. (14) repeatedly to reach convergence. 
F l ^° (t) = 1 for all i,j and t is a suitable starting con- 
dition, but the iteration itself can in practice be time- 
consuming and the calculation may not be tractable. 
Even if it is tractable, it almost certainly demands more 
effort than simply simulating the spread of an epidemic 
on the network of interest. There are some choices of the 
distributions r(r),s(T) for which the equations simplify 
and are more tractable — we examine two in Section V. 
Alternatively, one may be able to make useful approxi- 
mations in some cases. For instance, if /(r) is sharply 
peaked close to t = 0, as it is for many real diseases, 
then it may be reasonable to approximate F l< ~ J (t — r) 
in Eq. (14) by its value F^(t) at r = 0. Then (14) 
becomes 



F^{t) = 1 -p(t) + zp(t) ft f*-'(t), (22) 



where p(t) = J* /(r) dr. Hence the values of F 1 ^ 3 at 
different times decouple and the equations can be solved 
by a simple scalar iteration — no integrals need be per- 
formed. Although efficient, however, this approximation 
is usually only a good one in regions where F 1 ^ 3 (t) is rel- 
atively constant over the timescales typical of the disease 
progression represented by /(t), which means early and 
late times, but not in the crucial intermediate interval 
where most of the interesting behavior falls. 

Even in cases where the message passing equations are 
not a practical calculational tool, however, they can still 
be useful. In particular, they can tell us about the late- 
time limit of epidemics, including important quantities 
such as the total number of people infected by the dis- 
ease, and they allow us to calculate epidemic outcomes 
averaged over ensembles of networks such as the widely- 
studied configuration model. We look at these two ap- 
plications now in turn. 
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IV. LATE-TIME BEHAVIOR 



Taking the limit t — > oo in Eq. (14), we get 

/•OO 

F^'(oo) = l- / /(r) 1-z |~[ F j ^ l (oc) 
Jo ^ le/V(j)\i 



dr, (23) 



where we have assumed that /(r) is suitably small for 
large values of its argument. Writing F l4 ^° = F l ^ J (oo) 
for short and defining P = J /(t) dr, which is the total 
probability of transmission occurring between two ver- 
tices connected by an edge, we then find that 



= l-p + pz ] j F j< 

ieM( )\i 



(24) 



This again takes the form of a message passing calcula- 
tion, but now the messages passed are simple numbers, 
and hence the calculation can be performed quickly, even 
on networks that are not trees. Then the probability that 
a vertex is susceptible in the limit of late times satisfies 



P(Si)>z 11 F^i. 



(25) 



In the limit of late times there are no infected 
individuals — all have either recovered or never got sick 
in the first place — so P(Ri) = 1 — P(Si). Thus this cal- 
culation gives us an upper bound on the probability that 
any given individual ever contracts the disease or, if we 
sum over all vertices, an upper bound on the size of the 
disease outbreak. 

As has been discussed previously [24-27], the late-time 
limit of the SIR model is related to a correlated bond per- 
colation process on the corresponding directed transmis- 
sion network, the correlations arising because of variation 
in the time an individual takes to recover: if an individ- 
ual recovers quickly then the probability of transmission 
of the disease to any of its neighbors is small; if it takes 
a long time to recover the probability is correspondingly 
larger. Equations (24) and (25) can be considered to 
define a message passing algorithm for solving precisely 
this bond percolation problem on a general network. In 
this context, F l ^ J is a generating function in z for the 
number of vertices in the percolation cluster of vertex i 
that are reachable via vertex j, and P(Si) is a generating 
function for the overall sizes of the clusters. In recent un- 
published work, Shiraki and Kabashima [28] have given a 
message passing method for calculating percolation clus- 
ter sizes on trees and locally-tree-like networks, which is 
equivalent to the method reported here for the special 
case of a tree. 



V. EPIDEMICS ON CONFIGURATION MODEL 
NETWORKS 

Our method can also be used to calculate average prob- 
abilities of infection for ensemble models of networks. It 



is common in the study of processes on networks to look 
at not the behavior on a single network, but the aver- 
age behavior in an ensemble model defined as a proba- 
bility distribution over possible networks. The message 
passing formalism developed here allows one to calcu- 
late such average behaviors easily. We demonstrate this 
type of calculation using the configuration model, which 
is probably the most widely studied ensemble model of a 
network [29, 30]. 

In the configuration model one fixes the degree distri- 
bution of the network — meaning the fractions pk of ver- 
tices with each possible degree k — but in other respects 
connects vertices at random. Thus in calculating the be- 
havior of an epidemic on the configuration model there 
are two sources of randomness to average over. The first 
is the randomness in the dynamics of the disease, which 
is already built into our message passing formalism. The 
second is the randomness of the graph ensemble. 

Consider the average over the graph ensemble and con- 
sider an edge attached to vertex i. In different networks 
of the ensemble this edge will be attached to different 
vertices j at its other end and hence a different message 
H l ^ J will be transmitted down the edge. The ensem- 
ble average probability that vertex i has not yet been 
infected along the edge by time t is equal to the average 
of these messages over the set of networks. But, since 
every edge plays an identical role in the configuration 
model ensemble, the average message is the same for all 
edges i, j and hence we need calculate only one message 
to solve for the average behavior of the model. Let us 
denote this average message by Hi (t) . 

To calculate the average message, we need to aver- 
age Eq. (2) (or its equivalent, Eq. (14) for non-tree net- 
works), which requires us to average the product on the 
right-hand side of the equation. The average of such a 
product is not in general equal to the product of the av- 
erage message, which potentially makes the calculation 
more complicated. However, in the limit of large net- 
work size, configuration model networks have the crucial 
property of being locally tree-like, with the shortest cy- 
cles in the network being of length O(logn) and hence 
diverging as n — > oo. This means that the messages a 
vertex receives along each of its incident edges are inde- 
pendent in the large-n limit — in essence, we assume that 
correlations along a cycle of diverging length are irrele- 
vant in the large size limit. In this case, the average of 
the product of messages received by a vertex is equal to 
the product of the average. 

Averaging Eq. (2) over the ensemble and allowing for 
the fact that all messages are the same, the product 



n 



lGjV(j)\i 



H 3 in the equation now becomes simply a 



power [Hi(t)] k , where k is the so-called excess degree of j, 
i.e., its degree minus the edge between i and j, which has 
been removed because i is in the cavity state. The excess 
degree is distributed according to the excess degree dis- 
tribution q k ~ [k + l)pk+\/{k) [30] and, averaging over 
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this distribution, we find 

00 r ft 

1-/ /(^(l-Z^-T^dT 

fc=o L Jo 

= 1- /*/W[l-«G!i(iri(t-T))]dT, (26) 
Jo 

where G\(x) = ^ k qkX k is the probability generating 
function for the excess degree distribution. 

Similarly, from Eq. (3), the probability that a vertex 
of (ordinary) degree k is susceptible at time t is z[Hi(t)] k 
and the average probability of being susceptible is 

00 

P(S) = zY,Pk[Hi{t)] k = zG (Hi{t)), (27) 



fe=0 



where Go(x) = ^Z k PkX k is the generating function for 
the ordinary degree distribution p k . 

Again we can study the late-time behavior by letting 
t — > 00 and writing Hi = iZi(oo) to give 



and 



H 1 = l-p + pzG 1 (H 1 ), 



P(S) = zGo(ffi), 



(28) 



(29) 



where p = / °° /(r)dr. These two equations are precisely 
the standard equations for bond percolation on the con- 
figuration model [31] and highlight again the connection 
between the SIR model and percolation theory. The mes- 
sage Hi can be regarded as a generating function in z for 
the distribution of numbers of vertices reachable along 
an edge in bond percolation and P(S) is a generating 
function for the sizes of clusters. 



VI. EXAMPLES 

As a first example of the application of our formalism, 
consider what happens when the distributions r(r) and 
s(t) take the standard exponential form, corresponding 
to stochastically constant probabilities of infection with 
and recovery from disease. Specifically, we assume that 
s(t) = (3c~ i3t and r(r) = 7C~ 7T , where P and 7 are the 
rates of infection and recovery. Then /(r) = /3e~' /3+7 ^ T 
and, making the substitution t' = t — T, Eq. (26) becomes 

Hi{t) = 1 - /3c-^+^' /' c^+^'' [1 - zGi(Hi(t'))] dt'. 
Jo 

(30) 

Differentiating with respect to t, we then find that Hi 
satisfies 

^ = I3(f3 + 1 )e-^ t j t e^ t '[l - zGi(Hi(t'))] dt' 

- f3[l - zGi(Hi(t))] 
= 1 -{P + 1 )Hi{t)+PzGi{Hi{t)). (31) 



with the initial condition -ffi(O) = 1. This differential 
equation has the solution 



t 



Jl 



du 



7- {P + i)u + PzGi(u)' 



(32) 



And once we have H\(t) we can use Eq. (27) to calculate 
P(S) and subsequently P{I) and P(R). In Fig. 3 (top 
two frames) we show the form of the resulting solution 
for the particular choice of a Poisson degree distribution, 
along with the results of numerical simulations of epi- 
demics spreading on the same networks. As the figure 
shows, the analytic and numerical approaches agree well, 
and take the familiar form of an SIR outbreak with a 
brief peak in the number of infected individuals followed 
by a sharp decline and corresponding rise in the number 
of recovered individuals. 

But now consider a second choice that is quite differ- 
ent but perhaps more realistic. In this case we assume 
that individuals once infected do not become infectious 
immediately, passing through a latent period before de- 
veloping a transmissible infection, and also that infected 
individuals do not start recovering from disease imme- 
diately upon infection as in the exponential model, but 
remain infected for a certain length of time then recover. 
A simple choice displaying these two behaviors is the "top 
hat" function 



P 



■[6{T-T s )-6{T-T r )], 



(33) 



with r r > t s , where 9(t) is the Heaviside step function. 
In this expression r s is the time at which the infected 
individual becomes infectious, r r is the time at which 
they recover, and p, as before, is the total probability of 
transmission. 

Inserting this form into Eq. (26) and again differenti- 
ating gives 



dHi 
~dF 



P 



-Bit 



[6{t - T r )[l - zGi{Hi(t - T r ))] 

-T s )[l-zGi(Hi(t-T s ))}]. (34) 



where again -ffi(O) = 1. The lower two panels of Fig. 3 
show the solution of this equation for the same Poisson 
degree distribution as previously, and p, r r , and t s chosen 
so as to give the same mean time of transmission and 
total transmission probability as in the exponential case. 
Fixing the total transmission probability to be the same 
also fixes the long-time behavior to be the same, as can 
be seen in the figure. 

The two calculations exponential and "top hat" ver- 
sions of /(t) — nonetheless give quite different results. 
The epidemic peaks around the same time in each (about 
t = 6 in the plots), but more individuals are infected at 
any time in the exponential case and the epidemic lasts 
longer. Furthermore, the top hat case shows distinctive 
waves of infection, of period roughly equal to t s , sepa- 
rated by intervals of comparatively lower disease activ- 
ity. These waves are caused by the appearance of dis- 
tinct "generations" in the spread of the disease as the 
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FIG. 3: Fraction of the population infected (left) and recovered (right) as a function of time for two different choices of the 
parameters of the model. Calculations were performed on configuration model networks of n — 10 5 vertices and Poisson degree 
distribution with mean 3. In the top two panels infection and recovery times are exponentially distributed as described in the 
text, with /3 = | and 7 = §• In the bottom two panels f(r) takes the "top hat" form of Eq. (33), with r s = 0.8, r r = 1, 
and p — 0.8. The initial condition was z = 0.999 in each case. Solid lines in each panel are the predictions of the theory; circles 
are simulation results, averaged over 100 runs. 



first round of disease carriers passes infection to the sec- 
ond, who some time later pass it to the third, and so on. 
Such waves of infection are observed in many real-world 
diseases but are absent from models using a conventional 
exponential distribution of infection times (although they 
can be represented in a crude fashion by introducing ad- 
ditional disease states, as in the so-called SEIR model). 

For other choices of degree distribution, including 
power-law, uniform, and exponential distributions, the 
predictions are qualitatively similar by and large, and 
agree similarly well with simulation results. The shapes 
of the curves are, however, significantly altered by differ- 
ent choices of the parameters r r and t s in the top hat 
case: as the values of r r and r s become better separated 
the waves of infection become blurred and ultimately im- 
possible to distinguish. Conversely, the waves become 
more pronounced if r r and r s are chosen closer to one 
another. 



VII. CONCLUSIONS 

In this paper, we have studied the SIR model of epi- 
demic disease on a contact network, in a generalized 
form that allows for non-constant probabilities of infec- 
tion and recovery, by contrast with conventional SIR cal- 
culations. Abandoning constant probabilities obliges us 
also to abandon the traditional differential equation ap- 
proach to solving the model, but we have shown that the 
problem can be reformulated instead in the language of 
message passing. We have given a message passing cal- 
culation that is exact on networks that take the form of 
trees (or are locally tree-like, as in random graphs) and 
provides a rigorous bound on the probabilities of disease 
states on non-tree-like networks. 

We have demonstrated the application of our approach 
to the calculation of the late-time behavior of the gener- 
alized SIR model and to the calculation of average prop- 
erties of the model within the random graph ensemble 
known as the configuration model. One could in princi- 
ple extend the calculations to other random graph ensem- 
bles, such as random graphs with degree correlations [32] 
or random graphs with clustering [33], or to calculations 
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on single networks (i.e., not ensembles). 

The approach taken here can be applied to other dy- 
namical models on networks, such as the SI or SEIR mod- 
els, again yielding exact results on trees or tree-like net- 
works and rigorous bounds in the non-tree case, and it is 
possible the approach could also be applied to threshold 
models [34]. At the moment, it's unclear whether models 
such as the SIS model in which vertices can return to 
past states can be tackled in the message passing frame- 
work. The developments for the SIR model relied on our 
having an exact message passing solution on a tree. We 
have not yet been able to find a similar solution for the 
SIS model and so the development of a message passing 
method for this model remains an open problem. 



Because the functions fi are non-negative and monotonic 
in their first argument, the factors in brackets [. . .] are 
either both positive or both negative and hence the entire 
expression is non- negative for any x and y. Now let x and 
y be independent random variables, both with the same 
distribution as x\ and let us average (A. 3) over x and y. 
After rearranging we find that 

\ i=l ' xi \ i—2 ' xi 

The same argument can now be applied to the remaining 
functions / 2 , . . . , /„ in turn, to demonstrate that 
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n/<) >n<^ 



(A.5) 



= 1 ' X! i-\ 



and the equivalent result naturally holds for averages over 
any of the variables: 



Appendix: Chebyshev Integral Inequality 

Let fi(xi, . . . , Xk), • • • , fn{xi, . . . , Xk) be a set of n non- 
negative functions that are monotone decreasing or in- 
creasing in each of their k real- valued arguments for fixed 
values of the other arguments. (They can be increasing 
in one argument and decreasing in another.) Then it can 
be proved that 



[J/i(a;i,...,a; fc )) > ]J(Mx u . . . 7 x k )), (A.l) 



where the average is over any distribution of the inde- 
pendent variables x\, . . . , x k - The proof is as follows. 
Let (/) denote the partial average 



I 



f(xi,.. .,Xj,x j+ i, . . .,x k )P{xi) . ..P(xj) dx 1 . ..dxj, 

(A.2) 

which is a function of the remaining arguments Xj+i 
to Xk ■ Then consider the following product for arbitrary 
x and y 



[fi {x, x 2 , ■ ■ ■ , x n ) - /i (y, x 2 , ■ ■ ■ , x n )] 

■ n n 

Q fi{x, x 2 , . ■ . , x n ) - Y[ My, 



i=2 



(A.3) 



n/<> > n</<>,/ 

l—l ' X 3 l—l 



(A.6) 



The remainder of the proof proceeds by induction. As- 
sume that 



II ' ) II / ( A - 7 ) 

i— 1 ' xi...Xj i—1 



for j < k. Averaging both sides over one additional vari- 
able Xj + i gives 



> 



— ^ ' x±...Xj : Xj^.i * i=l 



n^),,..x 3 .) • (a.8) 



Xj + 1 



But (/i) ,...,(/„) is itself a set of monotone 

non- negative functions of the variables Xj+i, . . . , Xk- Ap- 
plying Eq. (A.6) to this set, we then find that 



i=l ' X-L...Xj + 1 i-\ 



Applying induction and using Eq. (A.5) as the base case, 
the result is now established. 
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